Analysis of monocyte data (SCZ vs CTR): most advanced model
Plotting the first QC results:
Density plot
Covariate correlation
#Covariate correlation
MEs_lin(new.covs,new.covs)#PC-covariate correlation (nomal p)
PCA_cov_cor_R(imputed.covs, filt.df)Next step: QC on vst-normalized expression values.
ComplexHeatmap::pheatmap(cor(filt.df,filt.df), color = heatmap.color.code)#Check influence of covariates on data variance after transformation
pca <- plotPCA.custom(as.matrix(trnsf.df), intgroup=c("status", "batch", "sex", "smoking"),
ntop = 50000, returnData=TRUE, metadata=imputed.covs, pc.1 = 1, pc.2 = 2)
PoV <- round(100 * attr(pca, "percentVar"))
PCAplot(pca, "batch", PoV.df=PoV, colors=c('orange','lightblue','lightgreen',"yellow"), pc.1 = 1, pc.2 = 2)
PCAplot(pca, "status", PoV.df=PoV, pc.1 = 1, pc.2 = 2)#Check influence of covariates on data variance after transformation
plotVarPart(vp)print("The model used is:")[1] “The model used is:”
dds.2@design~Mean.Per.Base.Cov. + Mapped + Exonic.Rate + Genes.Detected + rRNA.rate + RIN + BMI + age_clusters + batch + status
#PC-covariate correlation (nomal p)
PCA_cov_cor_R(imputed.covs, batch.rem)PCAplot(pca.cor, "batch", PoV.df=PoV.cor, colors=c('orange','lightblue','lightgreen',"yellow"), pc.1 = 1, pc.2 = 2)
PCAplot(pca.cor, "status", PoV.df=PoV.cor, pc.1 = 1, pc.2 = 2)We then go into the differential gene-expression analysis:
out of 24285 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 85, 0.35% LFC < 0 (down) : 111, 0.46% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 5) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results
As well as the number of differentially expressed genes at lfc </> -/+0.5 at padj < 0.1:
[1] 43
#DT::renderDT(data.frame(res.complex), "OH1.5A",scrollY=1000)
createDT(data.frame(res), "Status", scrollY=1000)EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
pCutoff = 0.1,
FCcutoff = 2,
labSize = 6,
xlim = c(-5,5),
ylim = c(-1,5),
legendPosition = 'right')#Plotting DEGs
ComplexHeatmap::pheatmap(batch.rem[rownames(batch.rem) %in% upreg.genes,], scale= "row",
cluster_rows = T, cluster_cols = T, annotation_legend = T, show_colnames = F, show_rownames = T,
legend = T, treeheight_row = 0, color = heatmap.color.code,
annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)ComplexHeatmap::pheatmap(batch.rem[rownames(batch.rem) %in% downreg.genes,], scale= "row",
cluster_rows = T, cluster_cols = T, annotation_legend = T, show_colnames = F, show_rownames = T,
legend = T, treeheight_row = 0, color = heatmap.color.code,
annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)#DT::renderDT(data.frame(res.complex), "OH1.5A",scrollY=1000)
paste("Reminder: we have", length(upreg.genes), "upregulated and", length(downreg.genes), "downregulated genes.")[1] “Reminder: we have 19 upregulated and 10 downregulated genes.”
createDT(DEG.enrich.res, "Enrichment", scrollY=1000)ComplexHeatmap::pheatmap(DEG.enrich, cluster_rows = F, cluster_cols = F, annotation_legend = F, show_colnames = T,
legend_breaks= c(0,-log(0.05), 5, 8, 10),
legend = T, color = colorRampPalette(RColorBrewer::brewer.pal(6,"RdPu"))(30),
legend_labels=c("0","2.995", "5", "8","-log(q) \n\n\n"))dotplot(GO.up, split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale="free")dotplot(GO.down, split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale="free")for (i in colnames(gene.panels_mrgd)){
print(i)
dat <- batch.rem[rownames(batch.rem) %in% na.omit(gene.panels_mrgd[,i]),]
plots <- PCA_cov_cor_R(imputed.covs, dat)
list.pca <- plotPCA.custom(as.matrix(dat), intgroup=c("status", "batch", "sex", "RIN"),
ntop = 50000, returnData=TRUE, metadata=imputed.covs, pc.1=1, pc.2=2)
list.pov <- round(100 * attr(list.pca, "percentVar"))
print(PCAplot(list.pca, "status", Shape="batch", PoV.df=list.pov, pc.1=1, pc.2=2))}[1] “Patir_microglia” [1] “Astrocytes” [1] “IFNy” [1] “NfKB” [1] “IFN.panel” [1] “LPS” [1] “IL4” [1] “GC.chronic” [1] “GC.acute.down” [1] “GC.acute.up”
#{r var.2, echo=TRUE, fig.height = 4, fig.width = 5.5, fig.align = "center"} #for (i in colnames(gene.panels_mrgd)){ # print(i) # dat <- batch.rem[rownames(batch.rem) %in% na.omit(gene.panels_mrgd[,i]),] # list.pca <- plotPCA.custom(as.matrix(dat), intgroup=c("status", "batch", "sex", "RIN"), # ntop = 50000, returnData=TRUE, metadata=imputed.covs, pc.1=1, pc.2=2) # list.pov <- round(100 * attr(list.pca, "percentVar")) # print(PCAplot(list.pca, "status", Shape="batch", PoV.df=list.pov, pc.1=1, pc.2=2))} #
for (i in colnames(gene.panels_mrgd)){
panel <- panel.list[[i]]
DE <- DE.list[[i]]
df <- panel[DE$log2FoldChange < -0.5 | DE$log2FoldChange > 0.5,]
df <- df[,order(as.numeric(colnames(df)), method="radix", decreasing=F)]
print(i)
ComplexHeatmap::pheatmap(as.matrix(df), scale= "row",
cluster_rows = T, cluster_cols = F, annotation_legend = T, show_colnames = F, show_rownames = T,
legend = T, treeheight_row = 0, color = heatmap.color.code,
annotation_col=annot_df, annotation_colors = annot_colors, fontsize_row = 8)
}[1] “Patir_microglia” [1] “Astrocytes” [1] “IFNy” [1] “NfKB” [1] “IFN.panel” [1] “LPS” [1] “IL4” [1] “GC.chronic” [1] “GC.acute.down” [1] “GC.acute.up”
(1A) Compound median lFC for each signature
#median lFC to point directions
ComplexHeatmap::pheatmap(-compoundlFC, color = heatmap.color.code)ComplexHeatmap::pheatmap(cbind(GEX.stat.cor[,1]), cluster_rows = T, cluster_cols = T, annotation_legend = F, show_colnames = T,
legend_breaks= c(0,0.01,0.05,0.1,0.1),
legend = T, color = colorRampPalette(RColorBrewer::brewer.pal(6,"RdPu"))(30),
legend_labels=c("0","0.01","0.05","0.1","Rho \n\n\n"),
display_numbers = as.matrix(matrix_pvalue[,1]))
ComplexHeatmap::pheatmap(cbind(GEX.stat.cor[,3]), cluster_rows = T, cluster_cols = T, annotation_legend = F, show_colnames = T,
legend_breaks= c(0,1,2.5,2.5),
legend = T, color = colorRampPalette(RColorBrewer::brewer.pal(6,"RdPu"))(30),
legend_labels=c("0","1","2.5","beta \n\n\n"),
display_numbers = as.matrix(matrix_pvalue[,2]))(2B) As differential PC expression
scatters.allclean.compoundScatters (3A) Gene program-Gene program interaction based on PCs
clean.compoundScatters2(3B) Data summarized as correlation on heatmap
ComplexHeatmap::pheatmap(rcorr(as.matrix(cbind(DEGscatter.df,"Status"=scatter.df[,1])), type="spearman")$r,
color = heatmap.color.code, fontsize_row = 8, fontsize_col = 8)ComplexHeatmap::pheatmap(z, cluster_rows = T, cluster_cols = T, annotation_legend = F, show_colnames = T,
legend_breaks= c(-2,-1,0,1,2,2),
legend = T, color = heatmap.color.code,
legend_labels=c("-2","-1","0","1","2","z-score \n\n\n"),
display_numbers = P)